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Abstract 



^ . The kinetic energy distribution function satisfying the Boltzmann equation 

1 is studied analytically and numerically for a system of inelastic hard spheres 



in the case of binary collisions. Analytically, this function is shown to have 



a similarity form in the simple cases of uniform or steady-state flows. This 
ON ■ determines the region of validity of hydrodynamic description. The latter is used 

1 to construct the phase diagram of granular systems, and discriminate between 

S clustering instability and inelastic collapse. The molecular dynamics results 

support analytical results, but also exhibit a novel fluctuational breakdown of 
O . mean-field descriptions. 

o 
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Granular media such as sand provide an attractive opportunity to revisit a number 
of topics in classical physics, and contribute new angles. In this paper we describe 
the granular phase diagram which we hope will be helpful to a broad community 
given the raising interest in granular systems pL The researchers interested in 
diluted granular gases, such as in astrophysical applications Q , and researchers who 
study, say, compaction of sand Q use different approaches. The phase diagram may 
represent a ground for communication. 

The phase of granular system depends on inelasticity of collisions, r, (restitution 
coefficient approximation, |Q, H |^]), particle density, p, particle size, a, system size, 

*J. Stat. Phys. (in press) 



L, and the observation time, t. The external influence of shaking, gravity, boundaries 
enters through the above parameters. Granular temperature, or any other charac- 
teristic of the rate of bulk-averaged motion doesn't appear on the phase diagram 
since there is no characteristic energy scale for hard-core interaction assumed here. 
As we shall see below, a combination, pLa^ 1 , where d is the system dimension, is 
playing the key role. It represents the average number of particles inside an imag- 
inary tube of length L and cross-section a d ^ 1 . Presently we think that the phase 
diagram consists of at least three regions, Fig. [|. In the region 1 — r <C (pLa^ 1 )^ 2 
the system is gas-like. In the region (pLa d ^ 1 )^ 2 « 1 - r « (pLa d ^ 1 )^ 1 the system 
is condensed but doesn't collapse, and for (pLa^ 1 )^ 1 <C 1 — r the system contains 
inelastically collapsed chains of particles, mixed with non-collapsing regions (see |7) 
for the description of collapse). Fig. is a snapshot from an two-dimensional event- 
driven Molecular Dynamics (MD) with 5000 particles in a circle the wall of which 
is maintained at a constant temperature. Particles undergo inelastic collisions with 
a constant coefficient of restitution (see Eqs. below). When a particle hits the 

wall it experiences diffusive angular scattering, while the distribution of the scattered 
velocity amplitude is Maxwellian, with the temperature being that of the wall. After 
a period of equilibration the system finds a state drawn in Fig. ^, which is in many 
respects a steady state. In Fig. || the gray scale linearly codes the relative number of 
collisions experienced by each particle during the previous 10 5 time steps, where black 
means high collision rate. One can clearly distinguish three different regions: (i) close 
to the wall there is a gas-like phase with low density where the mean free path is large, 
(ii) The region between the center of the bulk and the wall consists of closely packed 
particles. We consider this region as condensed phase. It is well separated from the 
gas like phase, (iii) Close to the center of the container we find a region which isn't 
distinguishable from the second region by "naked eye" . With the help of gray scale 
coding one observes collision chains, i.e. almost linearly arranged chains of particles 
which made major contribution to the previous 10 5 collisions, in comparison to their 
neighbors. Quantitatively, about 5 percent of the particles participated in about 96 
percent of the collisions. Sometimes we observed that very few (»5... 10) particles 



take almost all of the collisions in a certain time interval (not shown). The lengths 
of the appearing chains as well as their life times vary irregularly, their statistics will 
be discussed in detail elsewhere j8j. This figures serves as an illustration that in MD 
simulation one encounters all three regions of the phase diagram given in Fig. [I]. 

The role of observation time is not included here, and will be discussed separately 
below. Notwithstanding this classification, the condensed phases may be ordered in 
space (a crystal) or disordered (a glass). This different structure-based classification 
exists in the elastic limit. We now present the arguments used for constructing Fig. |l| 
beginning from the gas-like phase. 

As the parameter pLa^ 1 increases the first condensation or clustering transition 
of granular gas occurs. The above given criterion has been obtained by Goldhirsch 
and Zanetti || with the usage of granular hydrodynamics || . This description is 
based on the assumption of molecular chaos (^] and Maxwellian distribution functions 
for particle velocities Q which have been questioned and related to the concept of 
inelastic collapse ]Io[ |. In this paper we present the limits of validity of granular 
hydrodynamics, and then use the hydrodynamical clustering instability to argue that 
it occurs before the condition of inelastic collapse is satisfied. The latter is the line sep- 
arating collapsing and non-collapsing condensed phases, 1 — r ~ (pLa d ~ 1 )~ 1 , ^ [l0| |. 

Inelastic collapse, when a chain of particles experiences an infinite number of 
collisions in finite time as discussed by McNamara and Young @, is realized inside 
dense clusters, specified by granular hydrodynamics. Inside such condensed regions we 
have numerically observed collision chains. The relation of hydrodynamic instabilities 
of collision chains to their collapse and the importance of the upper line in Fig. ^ will 
be published elsewhere [||. 

We begin with diluted phase. The study of granular gases is greatly simplified by 
the fact that the binary collisions dominate and goes back to the works of Boltzmann 
and others JH], ^|. It is unclear a priori, whether inelasticity represents a regular 
or singular perturbation. The first part of this problem is the reduction of Boltz- 
mann equation to hydrodynamics, it can be studied independently from the known 
complications associated with hydrodynamical description, i.e. divergence of trans- 



port coefficients at high orders in spatial inhomogeneities and in low dimensions p"2| . 
The situation is somewhat analogous to kinetics of phonons where hydrodynamical 
reduction and second sound have a restricted range of applicability as compared to 
the kinetic equation More difficult problem is the validity of the Boltzmann 

equation itself, or the hypothesis of molecular chaos. 

The Boltzmann equation for inelastically colliding identical particles reads jy], |[ [j~4] 

dtf + vd s f = C(f,f) (1) 

where f(v,x,t) is the velocity distribution function, and C is the bilinear collision 

operator r r i ^ 

C(/>/) = J l^fif'-ffAlv-vildadvt . (2) 

The differential cross-section of two spheres of radius a is da = na 2 sin 2 6d6 / 2, where 
9 is the angle between the vectors q and v — v\ . An incoming collision event to the 
state {v, v\\ occurs between particles with velocities v ' and «/, 

1 + r 

v = v — q[q(v-V!)] (3) 

2r 

1 + r 

v[ = vi + — — q[q(v - Vi)] , (4) 

here q is a unit vector pointing from the center of the sphere 1 to the center of the 
sphere 2, and < r < 1 is the so-called coefficient of restitution, it models energy 

losses in the center of mass reference frame j?], |l5) . In the case r = 1 one recovers the 

usual elastic limit. 

We now consider a uniform cooling of granular gas without gravity inside an elastic 
two-dimensional circle when inelasticity is too small for developing clustering instabil- 
ity The problem is (hopefully!) isotropic and homogeneous and the distribution 
function, Jo, must depend only on the absolute value of particle velocity, v, i.e., on 
kinetic energy, E — mv 2 /2. When the initial distribution is forgotten one may search 
for a similarity solution 

= wAm) ' (5) 

where T(t) is a single scale of the kinetic energy. The normalization condition to 
the number density of particles n is / dz(f>(z) = 1. Eqs. (0, |) result in an integro- 



differential equation for function 4>(z) 



-Air) [ -(fr + z^J =^2 J dadztdag (zi) (6) 
x ^ —(^[(p' — cfxf>i ^ (z — Z\ — y/zzi cos a) 1 ^ 2 , 



and in an equation for T(t), 



dT/dt = -A(r)T 3/2 . (7) 



Here A(r) is the constant introduced by separation of variables, a is the angle between 
vectors v and vi, and g(z) is the density of states. Examination of the collision integral 

allows to get asymptotic form of <fi(z) at large z. Namely, it can be shown that the 

incoming term contains extra factor O (z -1 / 2 ) [this is contributed by the events 

(qv) oc z -1 / 2 ] and is small as compared to the outgoing term. The asymptotic form 

of Eq. (|) is -A(r)z{d<p/dz) - ^z 1 / 2 , up to a number. Therefore, 

ln#*)~V£M(r). (8) 

The moments of this function converge, and one can restrict oneself to the hydrody- 
namical description (Q) of the granular temperature which gives 



T{t) = T 



I + t -A{r)^/T 



(9) 



and T(t) oc t~ 2 at large t, Q. One can get this power law from dimensional ar- 
guments. A(r) remains undetermined, it can only be found from the full solution of 
Eq. (^). Asymptotically, A{r) oc (1 — r) when r is close to 1. The enhanced population 
of large energies, (||), dominates at z > l/^4 2 (r). For lower energies the distribution 
is Maxwellian. 

Two sets of MD results by an event-driven code are obtained with 5000 (20000) 
particles of unit radius in a circle with radius 130 (260). In both cases the surface frac- 
tion covered by particles is 50/169. Units of mass, length and time are arbitrary, initial 
velocity distribution is taken to be uniform in velocity in the square —1 < v x ,v v < 1. 
The employed model of collision is the same as given above, r = 0.999. Fig. || shows 
the temporal evolution of mean kinetic energy and the number of collisions occurred. 



The fit to the energy decay over the entire range is achieved with the help of Eq. ([)]). 
The distribution function is self-similar over four orders of energy decay as is seen in 
Figs. [| |. The latter fi gure shows the energy distribution function with the argument 
rescaled by (E(t)) = T for different times. One can see that the hydrodynamical 
description is quite precise. The system remains uniform and gas-like. 

However, at times t ~ 10 4 — 10 5 we observed appearance of a rarefied space in the 
center of the circle followed by "condensation" of granular gas on the wall. This is 
a novel type of transition, which is different from clustering || or collapse Q] . One 
can understand the growth of fluctuations using the following argument. The mean 
free path I is given by I = 1/na = a/c with the volume fraction c = na d in a d 
dimensional system with concentration n. The diffusivity of particles decreases in 
time like D(i) ~ (Sv)l ~ ly/T, where T is given in energy units. At large t (Eq. 

{tAf 



tl 



(10) 



Therefore, 

^2 



i(l-r) t(l 
and the diffusion length of particles in time t is 



I 2 (a/c) 
D~- = ) 1 > , (11) 



l D ^j Dmt ~^^{^), (12) 

where vq is the typical initial velocity. The number of particles inside the sphere of 
radius Id is N ~ nl d D = c (Id / a) d . These particles preserve a fluctuational collective 
velocity vq\J c(Id / a) d . Equating this velocity to thermal velocity Sv ~ y/T(i) ~ 
(a/c) / (t (1 — r)) yields the time when the collective velocity exceeds the fluctuational 
individual velocity 5 v 



d-3 
C 2 



^-TT^-^TTln^ 4 -±- , (13) 



vo(l- r )^ V 1 -^ 
given with logarithmic precision. Granular flow becomes fluctuationally supersonic 
at t 3> t c and molecular chaos assumption along with the Boltzmann equation and 
granular hydrodynamics are no longer applicable. The r.h.s. of Eq((l^) can be eval- 
uated, and the corresponding time is 5.0 x 10 4 . Given that the numerical prefactor 



in the estimate ([13]) isn't known, one finds a very satisfactory agreement between the 
estimate ( |i"3| ) and the time where hydrodynamical prediction deviates from the MD 
time-trace of kinetic energy density in Fig. |^. At time t c groups of particles of the 
size l c ~ [f c dtDlt)] 1 / 2 become special. Their mutual encounters lead to structures 
containing travelling shock waves. The difference between the systems with 5000 and 
20000 particles as seen in Fig. || is due to the fact that the regions of the size l c 
occupy a different fraction of the area, and in the second case it takes longer for 
inhomogeneities to evolve from the size l c up to the system size. The observed fluc- 
tuational transition is the reason why observation time may explicitly appear in the 
phase diagram. 

We now add a constant energy supply from the border, and allow the system 
to reach a steady state. The system is no longer uniform in space, and we assume 
that the corresponding length scales exceed the mean free path, I, everywhere. The 
similarity Ansatz Eq. (^|) with a non-uniform temperature T(x) can be used again. 
A study of the asymptotic form of the collision term shows that in the regions of 
the phase space where the detailed balance is absent the incoming term is again 
small in 0(z~ 1 ^ 2 ) times as compared to the outgoing term, and no global balance 
is possible. If 1 — r <C 1, the detailed balance approximately holds for the energies, 
z <C (1 — p) > above this range there are effectively no particles. Therefore, we find 
that the distribution function is almost Maxwellian for the energies less than T j (1 — r) 
and small otherwise, 

«w-{-*1 MIS. <"> 

All the moments converge, and one can justify the hydrodynamical reduction, which 
reads g | @ 

VP(p,T) = (15) 
Vq-B(r)T/r c (p,T) = , (16) 

where P is the granular pressure, and q is the thermal flux [Q, This system can 
be analyzed. To understand what happens as one increases inelasticity (or the number 
of particles) it is useful to consider the dilute limit approximation when q = X(T) VT, 



P = pT, where A(T) = 71T 1 / 2 is the thermal conductivity, B(r) oc (1 — r) is the 
energy fraction lost per collision. r c ~ 1 (n, T) = ^ipT 1 !" 1 is the time between collisions. 
We found that the hydrodynamical solution exists only below a threshold, 

a d-i B i/2 (r) f p dx < ^( 7l / 72 ) ln(L/a), (17) 



£d is a pure number. Our hydrodynamical analysis is valid when ([l7j ) is fulfilled; 
in the opposite case a particle condensate appears, and system has more than one 
phase. The left-hand side of (|l7|) is the generalization of the parameter pLa^ 1 for 
non-uniform case. Its region of applicability is more general than the derivation given 
above. Eq(|l7|) describes with logarithmic precision the same clustering transition as 
discussed by Goldhirsch and Zanetti ||. 

Returning to the question of validity of hydrodynamical reduction of the Boltz- 
mann equation we note that the complete system of equations for granular hydro- 
dynamics JTt| , ^, [l6| for density, linear and angular velocity density and granular 
temperature(s) has a shortcoming. It follows from Eqs. (PJT^) that in sufficiently in- 
elastic complex flows where changes occur in space and time the granular temperature 
cannot be introduced. Indeed, if we restrict ourselves with the first moment of the 
energy distribution in situations when this distribution changes its functional form, 
we wouldn't be able to close the hydrodynamical reduction. Therefore, in general, 
the full system of hydrodynamical equations [fl7j , cannot serve as a quantitative 
description. The situation is a bit easier, though, than in the kinetics of phonons: 
only kinetic coefficients are not exact, and the uncertainty is in numerical prefactors, 
which depend on local distribution function. To give an example about difficulties 
with phonons it is sufficient to recall that the heat transfer may be nonlocal Jl8[ . 

Sufficiently inelastic problems require a mixed kinetic-hydrodynamical description. 
In this case, similar to kinetics of phonons Jig] , one finds a reduced kinetic equation 
for the isotropic part, /o, of the distribution function 



dtfo-VJ[f ] = C(f ,f ), 



(18) 



where diffusion-like ffux J is defined as 

-l 

V/. (19) 



SC 



Sf 



1 

The inverse integral operator jj- is fixed by particle number constraint, i.e. the 



result of applying this operator to /o must have vanishing norm. Similar expressions 
arise for all dissipative coefficients. Despite the cumbersome appearance, Eq. ( |I§| ) 
together with the remaining hydrodynamic equations for density and momenta con- 
servations offers a reduction from nine to five-dimensional space of arguments. 

Our study of inelastic gas with binary collisions provided an opportunity to justify 
granular hydrodynamics for simple flows and for all flows in dilute and sufficiently 
elastic systems. This allows to discriminate between clustering instability discussed by 
Goldhirsch and Zanetti Q and inelastic collapse [Q. In complex inelastic flows the 
predictions of granular hydrodynamics are valid on the order of magnitude almost 
everywhere, and such precision is comparable to that in Fig. [I]. On the basis of 
this study we constructed the granular phase diagram and identified a novel type of 
supersonic fluctuational phase transition, leading to shock waves, which is different 
from clustering and collapse. 
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Figure 1: The granular phase diagram. See text for details. 




Figure 2: A snapshot from an event-driven Molecular Dynamics simulation of 5000 
hard spheres. The wall is kept at a fixed temperature. The gray scale codes the 
relative number of collisions experienced by particles during the previous 10 5 time 
steps. One distinguishes three regions: a gas-like state with low density and two high 
density regions, with and without collision chains (c.f. Fig. |l|). 




Figure 3: The kinetic energy E per particle over time t from MD simulations. Energy 
decay curves for systems with 5000 and 20000 particles are identical for a very long 
time (solid lines) The dashed line shows the prediction of Eq. [8] with T = 0.333 
and A = 0.01. The dash-dotted lines display the cumulative number of collisions per 
particle scaled by factors 10~ 8 (5000) and 4 x 10~ 8 (20000). 
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Figure 4: Evolution of the scaled energy distribution function in MD simulations. 
Different symbols used for different times as shown by insets. Statistically, the distri- 
bution is indistinguishable from Maxwcllian within the given range. 



